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ABSTRACT 

The  technique  of  windowing  has  been  often  used  in  the  implementation  of  the  waveform  relaj<- 
ations  for  solving  ODEs  or  time  dependent  PDEs.  Its  efficiency  depends  upon  problem  stiffness 
and  operator  spbtting.  Using  model  problems,  the  estimates  for  window  length  and  convergence 
rate  are  derived.  The  effectiveness  of  windowing  is  then  investigated  for  non-stiff  and  stiff  cases 
respectively.  It  concludes  that  for  the  former,  windowing  is  highly  recommended  when  a  large 
discrepancy  exists  between  the  convergence  rate  on  a  time  interval  and  the  ones  on  its  subintervals. 
For  the  latter,  windowing  does  not  provide  any  computational  advantage  if  mr  'line  features  are 
disregarded.  The  discussion  is  supported  by  experimental  results. 
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1.  Introduction.  The  waveform  relaxation  (WR)  method  v/as  originally  proposed 
for  solving  ordinary  differential  equations  arising  from  very  large  scale  integration 
(VLSI)  circuit  simulation  [5]  [9].  Unlike  conventional  timestepping  methods,  it  iter¬ 
atively  partitions  a  big  system  into  mutually  decoupled  subsystems,  and  then  solves 
each  subsystem  independently.  Different  discretizations  and  time  steps  are  allowed  for 
integrating  subsystems.  Based  on  its  nature,  the  method  has  been  proposed  as  a  mul¬ 
tirate  method  for  sequential  computing  or  a  parallel  method  on  advanced  computers 
[2]. 

Under  reasonable  assumptions  for  an  ODE,  the  WR  iteration  hcis  been  shown  to 
converge  superlinearly  on  finite  intervals  [1]  [5]  [8].  The  uniform  convergence  on  an 
interval  of  [0,  T]  is  reached  in  the  exponential  norm 

ll^lkr  :=  max  |e-^‘u(<)l,  {  >  0, 

which  implies  that,  for  many  problems,  the  WR  iteration  will  converge  much  faster 
on  short  intervals  than  on  longer  ones.  In  order  to  accelerate  the  convergence,  the 
technique  of  windowing  is  recommended,  in  which  the  interval  of  integration  is  split 
into  a  series  of  subintervals,  called  windows,  with  iteration  taking  place  on  each  window 
successively. 

The  length  of  window  is  of  practical  important  and  strongly  depends  upon  problem 
and  machine  involved.  The  general  guidance  for  its  selection  and  the  way  for  evaluating 
its  effectiveness  are  relatively  unknown  even  though  windowing  has  been  a  common 
practice  in  using  the  WR  method. 

The  estimates  for  time  windows  of  the  WR  iteration  were  studied  by  Leimkuhler 
and  Ruehli  for  RC  circuits  arising  as  simplified  models  of  a  VLSI  interconnect  [4].  Finer 
estimates  were  developed  by  Leimkuhler  for  a  model  linear  second-order  system  [3].  Us¬ 
ing  the  speed  of  splitting  and  weighted  spectral  radius  of  iteration  operator,  Leimkuhler 
estimated  the  abscissa  of  w-convergence,  which  then  provided  a  priori  estimate  for  the 
length  of  a  window  wherein  convergence  was  approximately  geometric  with  the  given 
rate  w.  His  approach  puts  emphasize  on  the  qualitative  comparisons  between  splittings. 

In  this  paper,  we  focus  on  the  time  dependency  of  the  approximation  error  and 
intimate  relation  between  the  WR  iteration  (or  dynamic  iteration  called  in  [6])  for 
time  dependent  problems  and  the  static  iteration  for  corresponding  steady-state  (or 
static)  problems.  For  certain  model  problems,  it  is  possible  to  separate  the  factor 
that  represents  early  sweeps  from  the  one  that  dominates  asymptotic  behavior.  Simple 
convergence  estimates  are  therefore  obtained.  The  estimates  and  the  results  observed 
in  the  experiments  are  compared  and  shown  to  have  good  agreements.  Based  on  these 
estimates,  the  effectiveness  of  windowing  is  discussed  for  non-stiff  and  stiff  problems 
respectively.  General  guidance  for  the  use  of  windowing  is  concluded  in  the  end. 

2.  Waveform  relaxation.  Using  a  first-order  linear  system 

du 

(1)  ~  +  Lu  =  f,  t>0,  u(0)  =  uo. 
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with  a  given  splitting  L  =  M  —  N,  the  WR  iteration  can  be  illustrated  by 


^  +  /,  t>0,  =  uo. 

at 

It  is  an  iterative  process  on  a  space  of  differential  functions.  The  functions  u^'’\  so 
called  “waveforms,”  are  then  discretized  for  numerical  integration.  The  continuous 
approximation  error  :=  —  u  satisfies 

(2)  t  >  0, 

where  5  is  a  linear  operator  on  (1  <  p  <  c»)  depending  upon  M  and  N, 

and  is  called  the  iUration  operator. 

Several  convergent  splittings  for  an  ODE  system  were  proposed  and  discussed  in 

[3]  and  [6].  We  shall  restrict  ourselves  to  splittings  that  resemble  Jacobi  splitting  and 
Gauss-Seidel  splitting  on  linear  systems  with  time  independent  coefficients,  as  described 
by  Eq.(l).  For  simplicity,  the  space  considered  in  this  work  is  C'‘*([0,  T],  C"),  the  space 
of  continuous  C"- valued  functions  in  [0,  T],  with  H-H  denotes  /“  norm  for  space  variables: 

IK0lh= 

and  III  ■  lllr  stands  for 


HMIIt  :=  m^lK<)||. 

The  notation  ||  *  ||  will  also  be  used  as  l°°  induced  matrix  and  operator  norm.  The  issue 
of  time  discretization  is  beyond  our  consideration,  for  which  the  reader  may  refer  to  [7]. 

3.  Convergence  estimates.  For  certain  type  of  problems  or  operator  splittings, 
the  factor  representing  early  phase  of  iterations  and  the  factor  dominating  asymptotic 
behavior  in  the  approximation  error  can  be  separated.  Laplace  transform  is  a  convenient 
tool  for  doing  this.  Applying  Laplace  transform  to  Eq.(2),  the  error  is  expressed  as 

eW(z)  =  5(z)e<’'-‘>(z)  =  5"(z)e<°)(z),  Rez  >  0, 


where 


S{z)  =  {zl -h  M)-^  N 

is  the  Laplace  transform  of  the  convolution  kernel  of  S.  Note,  5(0)  is  the  iteration 
operator  for  the  steady-state  problem  Lu  =  f  corresponding  to  Eq.(l). 

Theorem  1.  Let  L  be  split  as  L  =  M  —  N  with  M  =  dl,  d  >  0.  After  u  WR 
iterations,  the  error  is  bounded  by 

l|e'”>(<)ll  <  ».(*)  •  ll■S•(0)ll  •  ll|el“'|||T,  <  €  (OTl, 
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where 


=  g„(dt)  ■  l|S”(0)li  •  ll|e'”l|||T. 

Equation  (3)  can  be  easily  verified  by  induction.  □ 

It  is  interesting  to  examine  the  bounds  given  by  Theorem  1.  Note  that  (ll*5’'(0)|l)‘^’' 
is  time  independent  and  approximates  the  asymptotic  convergence  rate  either  for  ob¬ 
taining  the  steady-state  solution  or  solutions  over  long  time  intervals.  Function 
represents  the  time  dependency  of  the  error  and  dominates  the  convergence  behavior 
at  early  phaise  of  iteration  or  on  short  intervals.  It  is  a  monotone  increasing  function, 
bounded  by  0  <  gv{t)  <  1  with 


These  observations  on  the  error  bound  agrees  with  computational  experiences  that 
the  WR  iteration  converges  faster  on  short  intervals  than  on  longer  ones,  and  the 
convergence  rates  on  any  time  intervals,  including  infinite  interval,  are  at  least  as  good 
as  the  one  for  the  static  iteration. 

When  operator  L  has  constant  diagonal.  Theorem  1  actually  gives  an  error  bound 
for  Jacobi  WR  iteration.  The  next  two  theorems  will  give  results  for  Gauss-Seidel  WR 
iteration  on  model  problems. 

Theorem  2.  Let  L  =  M  —  N .  Assume  that  M  and  N  are  simultaneously  diago- 
nalizable  by  matrix  X,  and  all  eigenvalues  of  M  are  positive.  Then 

l|e<"'(()ll  <  9.W  ■  l|5''(0)ll  ■  amd{X)  ■  |||e'“>|||T,  t  6  10,7], 
where  d  is  the  largest  eigenvalue  of  M  and  cond{X)  =  ||A’“’||  •  11X11. 

Proof.  From  the  assumptions,  there  are  diagonal  matrices  Am  and  \n,  such  that 

M  =  X-^\mX,  N  =  X-^\sX, 

and 

Am  =  {A, (A/)},  Xi{M)  >  0  for  all  t. 

In  the  Laplace  domain, 

S"(z)  =  X-'{{zl  +  Am)-' Awrx  =  {X-'XlfA.1,X){X-'{zAJ^'  +  I)-X), 
leading  to 

■  /„(z/X,(Af)) 

eWM  =  5«(0)X-'  Xe«»)W. 

Applying  inverse  Laplace  transform,  one  obtains 

<  ll^'WII  •  ||X-'||  •  raaxj  /.(t)  ((  -  r /Ai(M))rfr| 

<  ||S”(0)||  ■  ||;f-||  •  raM(i,.(A,(M)()S|x,i|)  •  |||e<“>|||r 

;=1 

<,.(<«)'l|5''(0)l|-c<md(A-)-|||e"»|i|r.  □ 
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The  lexicographic  (or  forward  point)  Gauss-Seidel  WR  iteration  on  heat  equation 
with  periodic  boundary  condition  is  an  example  that  this  theorem  can  be  applied  (see 
[10]).  When  operator  L  has  the  form 


(4) 


L  =  d 


I  -B 
-R  I  ’ 


d>0, 


a  common  structure  when  using  a  red-black  ordering  on  certain  model  problems  in 
ODEs  or  in  time  dependent  PDEs,  the  error  in  Gauss-Seidel  WR  has  a  similar  upper 
bound. 


Theorem  3.  For  operator  L  represented  by  Eq.(4),  the  error  after  v  Gauss-Seidel 
WR  iterations,  expressed  as  e(*')(t)^  =  [eR^(<)^,CB^(<)^],  is  bounded  by 

(5)  l|eW(<)ll<to-.{<«)-|r(0)||-|||eg>|||T,  ielOTl- 


Proof.  For  Gauss-Seidel  splitting, 


M  =  d 


I  0 
-R  I 


,  and  N  =  d 


0  B 
0  0 


Laplace  transform  of  the  convolution  kernel  of  S  can  be  expressed  as 


S{z)  =  {zI  +  M)--^N  = 


0 


B 


[  ^  (*/4i)*  ^ 


P  =  RB,  Rez  >  0, 


which  yields 


S^{z)  = 


0  f2r,-i{zfd)BP'’-' 
[  0  f2v{zfd)P'’ 


and 


Back  to  the  time  domain. 


L  f,.{z/d)P"i‘'S\z) 


.(”)(()  = 


BP"-'  /„•"  -  rld)dr 


m, 


[  P"S^Mr)t'B\l-Tld)dT 


(0)/ 


Using  the  properties 


9vdp\{t)  <  9v{t)^  <>0,  u>0. 


max{l|BP-^||,ini}<|15’'(0)||, 
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and 


(note,  5"(0)  = 


O 

O  P'’ 


),  the  inequality  (5)  follows  immediately.  □ 


The  discussion  above  indicates  that  the  error  at  early  stage  of  the  WR  iteration  is 
controlled  by  function  g^.  Simple  convergence  estimates  are  therefore  derived  naturally 
from  this  function.  Given  a  convergence  rate  u>  and  an  iteration  number  v,  the  length 
of  a  window  of  u;- convergence  (see  [3])  can  be  estimated  by 

(6)  Tu,  ;=  mai{f  ;  {9mtdiv){dt)Y^'’  <  w}, 

where  mtd{-)  is  an  integer  function,  defined  by  the  operator  splitting.  For  instance, 
mtd{v)  =  V  and  mtd{v)  =  2v  —  \  for  the  cases  discussed  in  Theorem  1  and  Theorem  3 
respectively.  A  useful  variant  of  Eq.(6)  is  the  value 

(7)  :=  («„,d(.)(<ir))’7”, 

which  gives  an  estimated  average  convergence  rate  on  windows  of  length  T  in  first  v 
sweeps. 


Example  1.  Consider  the  ODE  system 

du  ^  ^ 

-  +  tu  =  0, 

where  L  =  (— 1,2,  —  1],  a  symmetric  tridiagonal  matrix  with  2  and  -1  on  main  and  off 
diagonal.  This  system  describes  the  nodal  voltage  of  a  linear  resistor-capacitor  (RC) 
network  [3]  [4].  Jacobi  WR  iteration  was  performed  with  randomly  generated  starting 
function  The  trapezoidal  rule  was  used  in  the  integration  with  conservative  time 
step  At  =  0.01  for  simulating  time-continuous  iteration.  For  a  given  number  of  sweeps 
V,  the  observed  convergence  rates  Uoh*  were  collected  as 


(max  max 
l<i<n  tg(0,rj 


l/v 


Figure  1  depicts  the  graph  of  for  5  Jacobi  iterations  (u  =  5)  together  with 
observed  data  marked  by  x’s.  Table  1  shows  more  detailed  comparison  between  the 
observed  convergence  rate  Uoba  and  estimated  rate  u>t  on  interval  [0,  T].  The  estimates 
TL’s  and  u^r’s  provided  by  Eq.(6)-(7)  are  surprisingly  close  to  the  observed  data.  Al¬ 
though  such  a  good  agreement  cannot  be  predicted  in  general,  reasonable  match  at 
early  phase  of  iteration  should  be  expected  if  similar  error  bounds,  «is  given  by  these 
theorems,  are  conjectured. 
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Table  1 

Observed  and  estimated  convergence  rate  on  [0,  T\ 


V 

Wo6*/wT 

T=0.2 

T=0.4 

T=0.6 

T=0.8 

T=1.0 

3 

.2271/.1994 

.3840/.3620 

.4920/.4939 

.5682/.6006 

.6227/. 6864 

5 

.1536/.1437 

.2793/.2691 

.3801/.3783 

.4622/.4730 

.5295/. 5550 

9 

.0000/.0927 

.1796/.1781 

.2555/.2568 

.3242/.3292 

.3852/.3956 

4.  Effectiveness  of  windowing.  The  theoretical  analysis  and  practical  experi¬ 
ments  have  revealed  that  the  effectiveness  of  the  WR  method  depends  highly  on  the 
stiffness  of  the  ODE  solved.  Not  surprisingly,  the  effectiveness  of  using  windows  in  the 
WR  iteration  is  also  closely  related  to  the  stiffness  of  the  system.  In  this  section,  we 
investigate  their  relation,  and  show  how  to  estimate  the  efficiency  of  windowing,  which 
then  results  the  general  guidance  for  the  implementation  of  the  WR  iteration. 

The  effectiveness  of  windowing  is  discussed  in  terms  of  computational  cost  or  op¬ 
eration  counts.  Following  concepts  and  notations  are  needed. 

Let  be  average  operation  counts  on  unit  windows  per  sweep.  For  a  given  error  tol¬ 
erance  e  >  0,  ut(c)  denotes  the  average  number  of  iterations  needed  for  the  convergence 
on  windows  of  length  T.  Using  Theorem  1,  vt(c)  can  be  estimated  as 

Ut(c)  :=  min{u  :  gmtd(v){dT)  <  c}. 

The  error  tolerance  e  will  be  dropped  whenever  the  context  is  clear.  The  total  com¬ 
putations  on  a  window  of  length  T  for  the  convergence  is  then  approximately  equal 
to 


C  :=  vtT fi. 

Note  that  Tfi  is  the  average  computations  on  a  window  of  length  T  per  sweep.  Let  this 
window  be  split  into  two  subwindows  of  length  Ti  and  T2,  Ti  -i-T-i  =  T,  with  the  WR 
iteration  taking  place  on  each  of  them  until  the  tolerance  level  is  reached  one  after  the 
other.  The  total  cost  then  satisfies 


vTiTiH  +  VTjTifi  <  vt'Th,  T'  =  max{Ti, r2}. 

Define  C-ujin  :=  vj'T g,.  C  and  (7w,n  will  be  referred  as  the  average  computations  on  an 
interval  of  length  T  without  and  with  windowing  respectively.  Since  vt'  <  vt,  we  have 

(8)  C^n  <  C. 


This  indicates  that  windowing  does  not  introduce  extra  computation. 

The  number  of  iterations  for  the  convergence  on  an  interval  of  length  T  is  closely 
related  to  the  stiffness  of  the  ODE  involved.  This  can  be  seen  from  the  behavior  of 
the  function  (see  Eq.(3)).  For  example,  an  approximation  for  the  case  stated  in 
Theorem  1  is 


g„{dT)  « 


{dTf-r 

vt\ 

7 


y 


leading  to 


That  is,  the  estimated  number  of  iterations  vt  increases  proportional  to  the  parameter 
d,  an  indicator  of  the  stiffness  of  the  ODE  solved. 

4.1.  Non-stifF  case.  Given  an  interval  of  length  T.  Following  above  arguments, 
the  number  of  the  WR  iterations  for  convergence  is  not  too  large  for  non-stiff  systems.  If 
the  interval  is  split  into  k  windows  of  equal  length  T'  =  T/k,  the  gain  or  the  percentage 
of  the  savings  of  using  windows  can  be  measured  by 

c,  _  ^  ~  _  vrTfi  -  kvT'T'fi 

~  C  ~  vtTh 

Vt  ~  vt' 

Vt 

If  VT  w  VT>,  Syjin  «  0,  not  much  computation  can  be  saved  by  windowing.  When 
Vt  >  VT',  using  windows,  vt  —  vt>  sweeps  of  WR  on  this  interval  of  length  T  are  likely 
to  be  reduced. 

Table  2  lists  the  experimental  results  on  Example  1  discussed  in  Section  3,  as  well 
as  corresponding  estimated  values.  The  entries  are  of  the  form  observed/estimated. 
The  numbers  listed  inside  parenthesis  in  column  2  are  average  number  of  iterations 
collected  on  subintervals  of  [0,  2]  with  length  T.  Windowing  was  used  on  intervals 
[0,r],  T  =  0.5,  1.0,  and  2.0  with  the  window  length  0.25.  The  iteration  was  terminated 
when  the  relative  error  was  in  the  order  of  c  =  l.e  —  7.  Again,  the  estimates  were 
in  good  agreement  with  the  observed  ones.  Windowing  reduced  the  computations  by 
25-60%.  It  clearly  suggests  that  windowing  is  quite  efficient  in  reducing  computations 
when  a  large  discrepancy  exists  between  the  convergence  rate  on  an  interval  and  the 
ones  on  its  subintervals. 


Table  2 

EfficUveness  of  windowing  (observed/estimated) 


Interval 

10,  r) 

No.  of  Iterations 
vt{  ave.  vt)/vt 

Computations 

^win 

% 

Without  Windowing 

With  Windowing 

[0,  0.25] 

8  (7)  /8 

2fi  /  2(i 

2fi  /  2n 

[0,  0.50] 

10  (9)  /lO 

bfi  /  bfi 

3.75/x  /  4/i 

25  /  20 

[0,  1.00] 

13  (12)  /13 

13/i  /  13^ 

7.25fi  /  S/i 

44  /  38 

[0,  2.00] 

17  (17)  /18 

34/i  /  36/i 

\Z.5fi  /  16|i 

60  /  56 

Remark.  Eq.(8)  and  the  observed  data  in  Table  2  seem  to  suggest  choosing  min¬ 
imum  window  length,  which  ironically  is  equivalent  to  the  step  size  used  in  the  time 
integration  of  subsystems.  In  this  situation,  the  WR  method  is  nothing  but  a  time¬ 
stepping  method  for  solving  ODEs.  However,  recall  that  the  method  is  proposed  as 
a  multirate  method  in  the  context  of  serial  computation  or  a  parallel  method  on  ad¬ 
vanced  computers.  For  the  former,  it  is  developed  for  problems  in  which  the  coupling 
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of  subsystems  is  relatively  loose  and  many  subsystems  allow  large  integration  steps. 
The  window  length  is  therefore  recommended  as  the  largest  step  size  used  in  the  time 
integration  of  subsystems.  For  the  later,  machine  characters  such  as  vector  length, 
communication  overhead,  play  important  roles.  The  study  in  this  paper  is  restricted  to 
the  mathematical  concerns  only. 

4.2.  Stiff  case.  In  this  situation,  the  WR  iteration  would  take  large  number  of 
iterations  to  converge  in  an  interval.  The  convergence  rate  very  likely  has  entered  the 
asymptotic  behavior.  Estimating  it  in  terms  of  function  alone  is  no  longer  adequate. 
From  the  error  bounds  given  by  Section  3,  the  rate  of  convergence  at  time  t  would  be 
dominated  by 

the  spectral  radius  of  5(0)  or  the  convergence  rate  of  the  related  static  iteration.  Since 
p(S(0))  is  time  independent,  the  convergence  rates  on  any  intervals  are  almost  identical, 
so  are  the  numbers  of  iterations  needed  for  the  convergence  on  those  intervals. 

Example  2.  Consider  the  heat  equation  on  the  unit  square  0  =  (0, 1)  x  (0, 1)  with 
Dirichlet  boundary  conditions 

ut  —  Au  =  0,  (t,  x)  €  (0,  T]  X  D 


u  =  0,  (<,x)  €  [0,T]  X  5ft,  u(0,i)  =  Uo(®)»  a;  €  ft. 

The  equation  was  first  discretized  in  space,  r^ulting  the  semi-discrete  problem 

(9)  ^  =  ®' 

where  is  the  five  point  difference  approximation  operator  to  the  Laplacian 


4  -1  , 

and  k  is  the  mesh  size  of  space  discretization.  The  red-black  Gauss-Seidel  WR  was  then 
implemented  on  the  system  (9)  over  time  interval  [0,7’].  The  iteration  was  terminated 
when  the  difference  between  the  uth  and  t;  —  1th  approximation 


reached  the  level  of  truncation  error  O(h^),  a  safe  stopping  criterion  proven  by  Nevan- 
linna  [8].  Table  3  shows  that  the  numbers  of  iterations  needed  on  different  time  intervals 
are  almost  the  same,  confirming  the  above  arguments.  As  is  discussed,  using  small  win¬ 
dows  for  this  problem  virtually  has  no  mathematical  advantage. 
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Table  3 

Number  of  iierutions 


WEM 

[0,0.125] 

[0,0.25] 

[0,0.5] 

[0,1.0] 

1/8 

16 

17 

17 

17 

1/16 

53 

60 

60 

60 

1/32 

208 

239 

239 

239 

5.  Conclusions.  In  this  paper,  the  convergence  estimates  such  as  window  length 
and  convergence  rate  are  developed  using  the  qualitative  comparison  between  the  WR 
iteration  and  the  corresponding  static  iteration.  The  effectiveness  of  windowing  tech¬ 
nique  for  the  WR  method  is  discussed.  The  results  proven  in  this  work  and  observed 
in  the  experiments  suggest  that,  for  uon-stiff  ODEs,  substantial  computations  can  be 
saved  by  windowing  when  a  large  discrepancy  exists  between  the  convergence  rate  on 
an  interval  and  the  ones  on  its  subinterval;  while  for  stiff  problems,  which  are  typically 
arisen  from  time  dependent  PDEs,  windowing  has  no  mathematical  advantage.  Thus 
for  stiff  systems,  the  selection  of  the  window  length  should  be  mainly  determined  by  the 
machine  features,  such  as  memory  capacity,  vector  length,  cache  size,  communication 
cost,  etc. 

Although  only  a  few  model  problems  are  considered  in  this  work,  similar  approach 
could  be  taken  for  some  generalized  problems.  The  guidance  concluded  above  certainly 
provide  helpful  information  in  the  implementation  of  the  WR  method  for  wide  class  of 
applications. 
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Fig  .  1 .  Observed  and  estimated  window  length 
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